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Abstract 

An iterative algorithm for the Helmholtz equation is presented. This 
scheme is based on the preconditioned conjugate gradient method for the normal 
equations. The preconditioning is one cycle of a multigrid method for the 
discrete Laplacian. The smoothing algorithm is red-black Gauss-Seidel and is 
constructed so it is a symmetric operator. The total number of iterations 
needed by the algorithm is independent of h. By varying the number of grids, 

O O 

the number of iterations depends only weakly on k when k h is constant. 
Comparisons with a SSOR preconditioner are presented. 
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Introduction 


We consider the Helmholtz equation 

Au + k^n^(x,y)u = f(x,y) (1.1) 

in a bounded two-dimensional region. In practice the region is unbounded in 
at least one direction. The unbounded region is then truncated with an 

artificial surface and a radiation boundary condition is imposed. In this 
paper we shall only consider Dirichlet and Neumann boundary conditions. 

Extensions to nonself-adjoint problems with local and nonlocal radiation 
boundary conditions verify the conclusions of this study. These results will 
be presented separately. The basic iterative method is presented in Bayliss, 
Goldstein, and Turkel [1]. There it is shown that a preconditioned conjugate 
gradient method is efficient for indefinite problem even when they are not 
self-adjoint. 

A discretization of (1.1) produces a large linear system of equations to 
be solved 

Au = b. (1.2) 

u is an approximation to the solution of (1.1) while b is determined by 
inhomogeneous data and boundary conditions. A is difficult to invert since 
it is not self-adjoint when radiation boundary conditions are used. In 
addition, the symmetric part of A is indefinite. Gaussian elimination is 
frequently used to solve (1.2). However, due to the large amount of storage 
needed, this method is limited to small k. As k increases, the solution 
becomes more oscillatory and more grid points are needed. In some simple 
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cases one requires k J h^ to be constant to achieve a fixed accuracy. In 
other cases the relationship of k to h can be worse than this [2]. 

To overcome this limitation an iterative method was proposed in [1]. The 
scheme is based on conjugate gradient for the normal equations. However, the 
normal equations have a condition number 0(h - ^) for fixed k, and so the 
iterations will converge slowly. We therefore precondition the matrix A by 
a partial inverse to the discrete Laplacian. Hence, we solve 

A'* A' x' = a'* b' (1.3) 

where 

A' = Q ^ A Q ^ , x' = x, b' = Q ^ b (1.4) 

and M~1 = Q -1 * Q“* is a partial inverse of Aq where A = Aq + k^ I. Thus, 
Aq is obtained by setting k = 0 in both the internal equations and the 
boundary conditions. We wish to choose M~^, so that the condition number 

"M 1 

of A' A ' is substantially reduced and M”* can be easily calculated. 

In [1], Q - ^ is constructed of one sweep of SOR. Hence, M - ^ is point 

SSOR applied to the Laplacian. In [1] it was also suggested that a more 
efficient preconditioning might be obtained from one cycle of multigrid. 
Since the multigrid algorithm requires only 0(n) operations, it is expected 
that the preconditioned problem has a convergence rate independent of h. A 
red-black Gauss-Seidel relaxation scheme is used since it is an efficient 
method for the Laplacian [4], Furthermore, it can lead to a symmetric 

preconditioner. 
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2. Iterative Algorithn 

The preconditioned conjugate gradient algorithm is described in detail in 

'k 

[1]. For this algorithm we need to compute A times a vector and A times 
a vector. In addition we need to evaluate acting on a vector where 

represents a multigrid cycle. We first briefly review the basic multigrid 
algorithm [3]. 

Given 

L U = f in fi, B U = g on 9ft , (2.1) 

we approximate (2.1) on the finest grid by 

t ® j -p\ m ..in m » / n o \ 

L U = F in G , BU=g on 9G , (2.2) 

The superscript in denotes an approximation on a grid with mesh h^ We also 
use auxiliary grids G* ,• • • with h^/hi = V2 • We denote by capital 
letters, the exact solution and by lower case letters, their approximation. 

1^ * denotes the weighting operator from to G^“^ and 1^ is the 

interpolation operator from to G^. 

The multigrid cycle begins with a "f ine-to-coarse" process which starts 
with one relaxation sweep on the finest grid. We then have u m which is an 
approximation to U m . We then transfer the residual 

R m = p m - L m u m (2.3) 


to the next coarser grid using 
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R m-1 = jitr-1 Rm> 
m 


( 2 . 4 ) 


On the grid 


G m_1 


we solve 


L ra-1 v m-l 


= R 


m-1 


( 2 . 5 ) 


where V m_1 


is an approximation to the fine grid error 


TT m TT m m 
V = U - u . 


( 2 . 6 ) 


This process is repeated to grids until we get to the coarsest grid. The 

size of this coarsest grid is arbitrary and will be varied in the result 
section. 

On the coarsest grid we use & relaxation sweeps rather than the one 

relaxation sweep used on the finer grids. To return to the finer grid we 

begin by using Si relaxation sweeps on the coarsest grid. The "coarse— to— 

2 

fine" process starts by calculating v new > which is a new approximation to 
V 2 , by 


2 


V 


new 


2 

V old 




( 2 . 7 ) 


2 

We then do one relaxation sweep to improve v new * This improved 
approximation is the correction to the next finer grid. This process is 
continued until we reach the finest grid, G m . This process consists of one 
cycle of multigrid. When used as a iteration process, this cycle is repeated 
until convergence is reached. Since we use multigrid only as a preconditioner 
we shall only use one cycle of multigrid. 
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In the basic scheme (see e.g., [5]) Gauss-Seidel relaxation is used. I k_1 

k 

is injection from to G^ ^ while is bilinear interpolation. For 

use as a preconditioning we require that the multigrid cycle be symmetric and 
positive definite, a characteristic which is not usually needed. Assuming a 
zero initial guess, the cycle is symmetric if all operations in the "fine-to- 
coarse" process are the adjoint of the operations in the "coarse-to-f ine" 
process. Specifically, 

(1) The fine-to-coarse relaxation must be the adjoint of the coarse-to- 
fine relaxation. 

(2) The coarse-to-f ine interpolation must be the adjoint of the fine-to- 
coarse residual weighting. 

(3) A symmetric operation must be performed on the coarsest grid. 

For the Laplace equation an efficient relaxation is Gauss-Seidel with red- 
black (RB) ordering [A]. It follows from property (1) that if RB is used when 
going to coarser grids, then BR must be used when returning to the finer 
grids. With this relaxation, full weighting in the fine— to— coarse residual 
transfer is preferable [5]. This agrees with property (2) for bilinear 
interpolation I^_^» The cycle is given schematically by Figure 1. 



Figure 1 
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Dirichlet boundary conditions present no difficulties. For the Neumann 
boundary conditions we introduce fictitious points and use full weighting [5]. 


3. Numerical Results 

The algorithm described above is valid for a general index of refraction 
n(x,y). In this section we shall only consider constant coefficients. 
Consider 

-Au - u = 0 0 _< x, y .< it (3.1) 

with boundary conditions 

u x (°, y ) = f(y) 

u (x,0) = 0 

y 

The analytic solution is given by 

u(x,y) = cos(^-)[ cos(k x (x - it )) - 2 cos(k x x)] , (3.3) 

9 9 

with k^ = k x + V4 • f(y), g(y) in (3.2) are given by those appropriate for 

(3.3). We use a uniform grid and (3.1) is approximated by 




0 . 


(3.4) 
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In Table I we present results using SSOR as the M” 1 preconditioner. We 
iterate until the residual is less than 10“ 6 . In Table II we present the 
results using one cycle of multigrid as the preconditioner. In this case we 
cycled down until we reached a 2x2 grid. We see that the number of cycles 
increases with finer mesh when SSOR preconditioning is used. Using multigrid, 
the number of iterations is independent of h, for fixed k. On the other 
hand, when h is fixed and k increases the number of iterations required by 
the multigrid preconditioner grows faster than that for the SSOR precondition- 
ing. For a realistic problem, k and h are not independent. The exact 
relationship between k and h to achieve a given accuracy depends on the 
geometry of the region, boundary conditions, and distance to eigenvalues [2]. 
A typical case requires k 3 h 2 to be constant for a fixed error tolerance. 
Hence, we wish to compare Tables I and II for this relationship. 

We first consider the effect of choosing the coarsest mesh to be finer 
than the 2x2 grid. In Table III we give the number of iterations for a given 
number of levels. Here one-level corresponds to 21 sweeps of Gauss-Seidel 
with red-black ordering without raultigrid and level = 2 is a two-level 
algorithm. In all cases, 4 relaxation sweeps RB followed by 4 BR, are used on 
the coarsest grid. We see that as k increases, it pays to use fewer levels. 
As we approach the Laplace equation, k = 0, we prefer to use all the possible 
levels. In other tests we varied the number of relaxation sweeps on the 
coarsest grid. As k increases, one should use fewer sweeps on the coarsest 
grid. We see from Table III that for k 3 h 2 fixed, the increase in iterations 
is a slowly increasing function of k. This growth is much slower than linear 
in k. However, when using a radiation condition, the growth seems to be 
closer to linear with k. Hence, this is much more efficient than SSOR. 
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The model problem considered is self-adjoint though indefinite. 
A. Bayliss (private communication) has done additional cases for a waveguide 
with a nonlocal radiation boundary condition which is no longer self-adjoint. 
Similar results were found which again show a slow growth with k when k 3 h 2 
is fixed. 


4. Conclusions 

A preconditioned conjugate gradient algorithm is presented for solving the 
Helmholtz equation. The preconditioning is based on one cycle of multigrid 
for the Laplace equation. The number of levels used is a function of k and 
h. As k increases or coarser meshes are used, few levels should be used. 
In this case one should also use fewer sweeps on the coarsest grid. 
Conversely, as k decreases or h gets smaller, one should use more levels 
and more relaxation sweeps on the coarsest grid. When the optimal choice of 
levels is made the number of iterations is a slowly (less than linear) growing 
function of k when k 3 h 2 is fixed. 



Table I. SSOR Preconditioning 


N 

8 

16 

32 

64 

k 





0.69 

27 

52 

108 

240 

1.39 

31 

61 

125 

272 

2.77 


69 

145 

333 

4.16 



166 

369 


Table II. 

Multigrid (to 2x2 grid) 

Preconditioning 


N 

8 

16 

32 

64 

k 





0.69 

9 

10 

10 

10 

1.39 

12 

14 

16 

16 

2.77 


33 

34 

31 

4.16 



80 

79 


Table III. 

Multigrid (to optimal, in computer time, 
coarsest grid) Preconditioning. 

In parenthesis is given the number of grids used 

N 

32 

64 

96 

k 




0.69 

10 (5) 

10 (6) 

11 (6) 

1.39 

12 (3) 

12 (4) 

15 (4) 

2.77 

15 (2) 

17 (3) 

19 (4) 

4.16 

22 (1) 

24 (2) 

21 (3) 

5.56 

21 (1) 

22 (2) 

29 (3) 
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